#define CATCH_CONFIG_MAIN
#include "catch.hpp" 
#include "discreteOscillators/discreteOscillators.h"
#include <iostream>

auto equal = [](auto x, auto y, double tol = 1e-5){return x == Approx(y).epsilon(tol);};

TEST_CASE( "Discrete oscillator treatment" ){
  GIVEN( "Test material" ){
    double lambda_s = 0.26460498561058793, temp = 296.0,
    twt = 0.0, tbeta = 0.5, effectiveTemp = 572.61028482016525;

    WHEN( "Single oscillator is present" ){
      std::vector<double> oscEnergies {0.3}, oscWeights {0.5};
      auto oscEnergiesWeights = ranges::view::zip(oscEnergies,oscWeights);
      AND_WHEN( "Few alpha, beta values spanning 0.1 -> 10's" ){
        std::vector<double> 
          alpha = {0.1, 1.0, 10},
          beta  = {0.0, 5.0, 50},
          sab { 6.28055E-3, 1.48422E-3, 1.78104E-29, 5.10424E-2, 1.52272E-2, 
          4.7449E-20, 7.01340E-2, 8.69301E-2, 8.7082E-10}, 
          correct { 6.25390E-3, 1.47793E-3, 0.00000000, 4.89179E-2, 1.45936E-2, 
          3.1085E-11, 4.58457E-2, 5.68386E-2, 8.52874E-5 };
          discreteOscillators( lambda_s, twt, tbeta, alpha, beta, temp, 
                  oscEnergiesWeights, effectiveTemp, sab );
        THEN( "scattering law matrix is correctly changed" ){
           REQUIRE(ranges::equal(sab,correct,equal));
        } // THEN
      } // AND WHEN
      AND_WHEN( "Very small alpha, beta values" ){
        std::vector<double> 
        alpha = {1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2},
        beta  = {1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2},
        sab { 6.42751E-7, 6.42911E-7, 6.42780E-7, 6.42911E-7, 6.43074E-7, 
        6.44381E-7, 6.46015E-7, 6.59084E-7, 3.21339E-5, 3.21419E-5, 3.21354E-5, 
        3.21419E-5, 3.21501E-5, 3.22154E-5, 3.22971E-5, 3.29504E-5, 6.42738E-6, 
        6.42898E-6, 6.42767E-6, 6.42898E-6, 6.43061E-6, 6.44368E-6, 6.46001E-6, 
        6.59070E-6, 3.21339E-5, 3.21419E-5, 3.21354E-5, 3.21419E-5, 3.21501E-5, 
        3.22154E-5, 3.22971E-5, 3.29504E-5, 6.42604E-5, 6.42764E-5, 6.42633E-5, 
        6.42764E-5, 6.42927E-5, 6.44234E-5, 6.45867E-5, 6.58933E-5, 3.21004E-4, 
        3.21084E-4, 3.21019E-4, 3.21084E-4, 3.21166E-4, 3.21819E-4, 3.22634E-4, 
        3.29161E-4, 6.41267E-4, 6.41427E-4, 6.41296E-4, 6.41427E-4, 6.41590E-4, 
        6.42893E-4, 6.44523E-4, 6.57560E-4, 3.17680E-3, 3.17760E-3, 3.17695E-3, 
        3.17760E-3, 3.17840E-3, 3.18486E-3, 3.19292E-3, 3.25746E-3 },
        correct { 6.42751E-7, 6.42911E-7, 6.42780E-7, 6.42911E-7, 6.43074E-7, 
        6.44381E-7, 6.46015E-7, 6.59084E-7, 3.21332E-5, 3.21412E-5, 3.21347E-5, 
        3.21412E-5, 3.21494E-5, 3.22147E-5, 3.22964E-5, 3.29497E-5, 6.42735E-6, 
        6.42895E-6, 6.42764E-6, 6.42895E-6, 6.43058E-6, 6.44365E-6, 6.45999E-6, 
        6.59067E-6, 3.21332E-5, 3.21412E-5, 3.21347E-5, 3.21412E-5, 3.21494E-5, 
        3.22147E-5, 3.22964E-5, 3.29497E-5, 6.42576E-5, 6.42736E-5, 6.42606E-5, 
        6.42736E-5, 6.42900E-5, 6.44206E-5, 6.45839E-5, 6.58905E-5, 3.20936E-4, 
        3.21016E-4, 3.20951E-4, 3.21016E-4, 3.21098E-4, 3.21750E-4, 3.22566E-4, 
        3.29091E-4, 6.40995E-4, 6.41154E-4, 6.41024E-4, 6.41154E-4, 6.41317E-4, 
        6.42620E-4, 6.44249E-4, 6.57280E-4, 3.17006E-3, 3.17085E-3, 3.17020E-3, 
        3.17085E-3, 3.17165E-3, 3.17809E-3, 3.18614E-3, 3.25055E-3 };
        discreteOscillators( lambda_s, twt, tbeta, alpha, beta, temp, 
                oscEnergiesWeights, effectiveTemp, sab );
        THEN( "scattering law matrix is correctly changed" ){
          REQUIRE(ranges::equal(sab,correct,equal));
        } // THEN
      } // AND WHEN
    } // WHEN

    WHEN( "Two oscillators present" ){
      std::vector<double> oscEnergies {0.02, 0.05}, oscWeights {0.2, 0.3};
      auto oscEnergiesWeights = ranges::view::zip(oscEnergies,oscWeights);
      AND_WHEN( "Few alpha, beta values spanning 0.1 -> 10's" ){
        std::vector<double> 
          alpha = {0.1, 1.0, 10},
          beta  = {0.0, 5.0, 50.0},
          sab { 6.28055E-3, 1.48422E-3, 1.78104E-29, 5.10424E-2, 1.52272E-2, 
                4.74495E-20, 7.01340E-2, 8.69301E-2, 8.7082E-10 },
          correct  { 1.43400E-2, 1.49937E-3, 0.00000E+00, 6.46197E-2,
                     1.71145E-2, 1.05458E-19, 9.47556E-3, 4.65210E-2,
                     2.20689E-8 };
          discreteOscillators( lambda_s, twt, tbeta, alpha, beta, temp, 
                  oscEnergiesWeights, effectiveTemp, sab );
        THEN( "scattering law matrix is correctly changed" ){
          REQUIRE(ranges::equal(sab,correct,equal));
          } // THEN
      } // AND WHEN
      AND_WHEN( "Very small alpha, beta values" ){
        std::vector<double> 
        alpha = {1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2},
        beta  = {1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2},
        sab { 6.42751E-7, 6.42911E-7, 6.42780E-7, 6.42911E-7, 6.43074E-7, 
        6.44381E-7, 6.46015E-7, 6.59084E-7, 3.21339E-5, 3.21419E-5, 3.21354E-5, 
        3.21419E-5, 3.21501E-5, 3.22154E-5, 3.22971E-5, 3.29504E-5, 6.42738E-6, 
        6.42898E-6, 6.42767E-6, 6.42898E-6, 6.43061E-6, 6.44368E-6, 6.46001E-6, 
        6.59070E-6, 3.21339E-5, 3.21419E-5, 3.21354E-5, 3.21419E-5, 3.21501E-5, 
        3.22154E-5, 3.22971E-5, 3.29504E-5, 6.42604E-5, 6.42764E-5, 6.42633E-5, 
        6.42764E-5, 6.42927E-5, 6.44234E-5, 6.45867E-5, 6.58933E-5, 3.21004E-4, 
        3.21084E-4, 3.21019E-4, 3.21084E-4, 3.21166E-4, 3.21819E-4, 3.22634E-4, 
        3.29161E-4, 6.41267E-4, 6.41427E-4, 6.41296E-4, 6.41427E-4, 6.41590E-4, 
        6.42893E-4, 6.44523E-4, 6.57560E-4, 3.17680E-3, 3.17760E-3, 3.17695E-3, 
        3.17760E-3, 3.17840E-3, 3.18486E-3, 3.19292E-3, 3.25746E-3 },
        correct  { 6.42745E-7, 6.42905E-7, 6.42775E-7, 6.42905E-7, 6.43069E-7, 
        6.44376E-7, 6.46009E-7, 6.59078E-7, 3.21196E-5, 3.21276E-5, 3.21211E-5, 
        3.21276E-5, 3.21358E-5, 3.22011E-5, 3.22827E-5, 3.29358E-5, 6.42681E-6, 
        6.42841E-6, 6.42710E-6, 6.42841E-6, 6.43004E-6, 6.44311E-6, 6.45944E-6, 
        6.59012E-6, 3.21196E-5, 3.21276E-5, 3.21211E-5, 3.21276E-5, 3.21358E-5, 
        3.22011E-5, 3.22827E-5, 3.29358E-5, 6.42034E-5, 6.42194E-5, 6.42064E-5, 
        6.42194E-5, 6.42357E-5, 6.43663E-5, 6.45294E-5, 6.58349E-5, 3.19784E-4, 
        3.19863E-4, 3.19798E-4, 3.19863E-4, 3.19945E-4, 3.20597E-4, 3.21415E-4, 
        3.28138E-4, 7.14086E-4, 7.14265E-4, 7.14119E-4, 7.14265E-4, 7.14451E-4, 
        7.16049E-4, 7.18329E-4, 7.48380E-4, 7.16740E-3, 7.16917E-3, 7.16773E-3, 
        7.16917E-3, 7.17098E-3, 7.18566E-3, 7.20446E-3, 7.37363E-3 };
        discreteOscillators( lambda_s, twt, tbeta, alpha, beta, temp, 
                oscEnergiesWeights, effectiveTemp, sab );
        THEN( "scattering law matrix is correctly changed" ){
          REQUIRE(ranges::equal(sab,correct,equal));
        } // THEN
      } // AND WHEN
      AND_WHEN( "Rather large alpha, beta values" ){
        std::vector<double> 
        alpha1 = {60, 100, 120},
        beta1  = {0,  100, 120},
        sab1 { 1.70896E-4, 5.21119E-8, 7.0407E-11, 2.39298E-6, 1.37913E-4, 
        2.23740E-6, 2.95504E-7, 1.14848E-3, 4.65867E-5 },
        correct1 { 1.36280E-8, 1.17072E-7, 2.12854E-8, 1.1482E-13, 2.15253E-6, 
        4.70458E-6, 8.8005E-16, 1.16734E-6, 5.99689E-6 },

        alpha2 {60, 100, 120, 150},
        beta2  {0,  100, 120, 150},
        sab2 { 1.70896E-4, 5.21119E-8, 7.0407E-11, 8.5425E-16, 2.39298E-6, 
        1.37913E-4, 2.23740E-6, 1.04788E-9, 2.95504E-7, 1.14848E-3, 4.65867E-5, 
        8.40991E-8, 1.31808E-8, 7.62407E-3, 9.66387E-4, 9.37698E-6},
        correct2 { 1.36280E-8, 1.17072E-7, 2.12854E-8, 1.1333E-11, 1.1482E-13, 
        2.15253E-6, 4.70458E-6, 4.31086E-7, 8.8005E-16, 1.16734E-6, 5.99689E-6, 
        4.54521E-6, 8.7052E-19, 1.97621E-7, 2.80726E-6, 2.03141E-5},

        alpha3 {10,30,60,100,120,150},
        beta3  {0, 50,100,120,150},
        sab3 { 7.01340E-2, 8.7082E-10, 1.0288E-23, 2.3337E-30, 0.00000000, 
        5.38240E-3, 4.48589E-5, 8.8775E-14, 5.0053E-18, 5.6793E-25, 1.70896E-4, 
        6.77023E-3, 5.21119E-8, 7.0407E-11, 8.5425E-16, 2.39298E-6, 2.85699E-2, 
        1.37913E-4, 2.23740E-6, 1.04788E-9, 2.95504E-7, 2.26650E-2, 1.14848E-3, 
        4.65867E-5, 8.40991E-8, 1.31808E-8, 8.35859E-3, 7.62407E-3, 9.66387E-4, 
        9.37698E-6},
        correct3 { 6.17897E-3, 1.34999E-8, 0.00000000, 0.000000000, 0.00000000, 
        2.36273E-5, 2.10282E-4, 3.7982E-10, 1.2303E-13, 0.00000000, 1.74719E-8, 
        4.57293E-4, 8.64524E-6, 7.73706E-8, 1.1333E-11, 1.8790E-13, 8.18326E-6, 
        3.00123E-4, 3.55855E-5, 4.32883E-7, 1.0024E-15, 3.77351E-7, 2.03275E-4, 
        6.32648E-5, 4.70616E-6, 8.7052E-19, 2.06917E-9, 4.03784E-5, 4.21806E-5, 
        2.22134E-5};

        discreteOscillators( lambda_s, twt, tbeta, alpha1, beta1, temp, 
                oscEnergiesWeights, effectiveTemp, sab1 );
        discreteOscillators( lambda_s, twt, tbeta, alpha2, beta2, temp, 
                oscEnergiesWeights, effectiveTemp, sab2 );
        discreteOscillators( lambda_s, twt, tbeta, alpha3, beta3, temp, 
                oscEnergiesWeights, effectiveTemp, sab3 );

        THEN( "scattering law matrix is correctly changed" ){
          REQUIRE(ranges::equal(sab1,correct1,equal));
          REQUIRE(ranges::equal(sab2,correct2,equal));
          REQUIRE(ranges::equal(sab3,correct3,equal));
        } // THEN
      } // AND WHEN
    } // WHEN

    WHEN( "Hefty number of oscillators present" ){
      std::vector<double> oscEnergies {0.01, 0.02, 0.05, 0.10, 0.20, 0.50}, 
                          oscWeights  {0.10, 0.10, 0.10, 0.10, 0.05, 0.05};
      auto oscEnergiesWeights = ranges::view::zip(oscEnergies,oscWeights);
      std::vector<double> 
      alpha1 {0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0},
      beta1  {0.0, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0},
      sab1 { 6.28055E-3, 3.47479E-3, 3.20577E-3, 3.99027E-3, 5.46999E-3, 
      1.48422E-3, 3.71333E-6, 7.25300E-9, 1.22741E-2, 6.81092E-3, 6.29090E-3, 
      7.82932E-3, 1.07336E-2, 2.98207E-3, 1.50949E-5, 5.91423E-8, 2.86328E-2, 
      1.60297E-2, 1.48572E-2, 1.84831E-2, 2.53465E-2, 7.53727E-3, 9.85209E-5, 
      9.73723E-7, 5.10424E-2, 2.89998E-2, 2.70325E-2, 3.36090E-2, 4.61129E-2, 
      1.52272E-2, 4.17548E-4, 8.37966E-6, 8.12099E-2, 4.75144E-2, 4.47911E-2, 
      5.56259E-2, 7.64146E-2, 3.03282E-2, 1.80045E-3, 7.44937E-5, 1.03383E-1, 
      6.59642E-2, 6.42001E-2, 7.95426E-2, 1.09831E-1, 6.62828E-2, 1.16676E-2, 
      1.31603E-3, 7.01340E-2, 5.12789E-2, 5.22934E-2, 6.47152E-2, 9.04794E-2, 
      8.69301E-2, 3.70511E-2, 9.46606E-3, 3.80120E-2, 3.13168E-2, 3.31684E-2, 
      4.11020E-2, 5.83713E-2, 7.69054E-2, 5.64106E-2, 2.39617E-2},
      correct1 { 5.67654E-3, 3.37158E-3, 1.64208E-1, 3.94132E-3, 5.31866E-3, 
      1.49950E-3, 8.46717E-6, 9.75667E-8, 1.01370E-2, 6.42619E-3, 2.71190E-1, 
      7.64226E-3, 1.01732E-2, 3.05814E-3, 3.57861E-5, 5.95960E-7, 1.90860E-2, 
      1.40460E-2, 3.93691E-1, 1.74244E-2, 2.25182E-2, 8.10371E-3, 1.78946E-4, 
      1.97549E-6, 2.71109E-2, 2.29336E-2, 3.51862E-1, 2.97806E-2, 3.77145E-2, 
      1.71846E-2, 7.50854E-4, 1.55913E-5, 3.29357E-2, 3.08740E-2, 1.95372E-1, 
      4.21364E-2, 5.44259E-2, 3.45691E-2, 3.07956E-3, 1.37933E-4, 2.54173E-2, 
      2.57397E-2, 5.04913E-2, 3.65298E-2, 5.14296E-2, 5.86794E-2, 1.67346E-2, 
      2.32145E-3, 6.99858E-3, 7.40993E-3, 8.79803E-3, 1.09541E-2, 1.69986E-2, 
      2.89090E-2, 2.06173E-2, 7.00193E-3, 1.48562E-3, 1.62702E-3, 1.95033E-3, 
      2.52931E-3, 4.18802E-3, 1.04260E-2, 1.48545E-2, 9.61884E-3 },

      alpha2  {0.1, 0.5, 1.0, 5.0, 10.0, 15.0},
      beta2   {0.0, 0.5, 1.0, 5.0, 10.0, 15.0},
      sab2    {6.28055E-3, 3.20577E-3, 3.99027E-3, 1.48422E-3, 3.71333E-6, 
      7.25300E-9, 2.86328E-2, 1.48572E-2, 1.84831E-2, 7.53727E-3, 9.85209E-5, 
      9.73723E-7, 5.10424E-2, 2.70325E-2, 3.36090E-2, 1.52272E-2, 4.17548E-4, 
      8.37966E-6, 1.03383E-1, 6.42001E-2, 7.95426E-2, 6.62828E-2, 1.16676E-2, 
      1.31603E-3, 7.01340E-2, 5.22934E-2, 6.47152E-2, 8.69301E-2, 3.70511E-2, 
      9.46606E-3, 3.80120E-2, 3.31684E-2, 4.11020E-2, 7.69054E-2, 5.64106E-2, 
      2.39617E-2},
      correct2 { 5.70932E-3, 1.32089E-1, 3.90917E-3, 1.47667E-3,7.18819E-6, 
      9.59641E-8, 1.94383E-2, 3.18818E-1, 1.70866E-2, 7.60852E-3, 1.77424E-4, 
      1.97549E-6, 2.75519E-2, 2.88058E-1, 2.92662E-2, 1.55401E-2, 7.50322E-4, 
      1.55909E-5, 2.42915E-2, 4.57593E-2, 3.57510E-2, 5.07956E-2, 1.66875E-2, 
      2.30261E-3, 6.64775E-3, 8.49761E-3, 1.09801E-2, 2.60825E-2, 2.04368E-2, 
      7.00483E-3, 1.40950E-3, 1.87368E-3, 2.45477E-3, 9.71456E-3, 1.46367E-2, 
      9.61785E-3},

      alpha3 {5.0, 10.0, 15.0},
      beta3 {5.0, 10.0, 15.0},
      sab3 { 6.62828E-2, 1.16676E-2, 1.31603E-3, 8.69301E-2, 3.70511E-2, 
      9.46606E-3, 7.69054E-2, 5.64106E-2, 2.39617E-2 },
      correct3 { 4.83392E-2, 1.66713E-2, 2.26650E-3, 3.18367E-2, 2.06018E-2, 
      7.00440E-3, 1.71694E-2, 1.53809E-2, 9.62749E-3};

      discreteOscillators( lambda_s, twt, tbeta, alpha1, beta1, temp, 
                oscEnergiesWeights, effectiveTemp, sab1 );
      discreteOscillators( lambda_s, twt, tbeta, alpha2, beta2, temp, 
                oscEnergiesWeights, effectiveTemp, sab2 );
      discreteOscillators( lambda_s, twt, tbeta, alpha3, beta3, temp, 
                oscEnergiesWeights, effectiveTemp, sab3 );


      THEN( "scattering law matrix is correctly changed" ){
        REQUIRE(ranges::equal(sab1,correct1,equal));
        REQUIRE(ranges::equal(sab2,correct2,equal));
        REQUIRE(ranges::equal(sab3,correct3,equal));
      } // THEN
    } // WHEN
  } // GIVEN

  GIVEN( "Simple water example" ){
    double temp = 400.0, lambda_s = 0.38283816, effectiveTemp = 600.970054,
    twt = 0.055556, tbeta = 0.444444;
    std::vector<double> oscEnergies {0.205,0.48}, oscWeights {0.166667,0.333333};
    auto oscEnergiesWeights = ranges::view::zip(oscEnergies,oscWeights);

    WHEN( "larger alpha, beta grid" ){
      std::vector<double> 
      alpha{0.01, 0.02, 0.04, 0.06, 0.08, 0.20, 0.40, 0.60, 0.80, 2.00, 4.00, 
      6.00, 8.00, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0 },
      beta {0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 
      0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 2.00, 4.00, 6.00, 
      8.00, 10.0, 15.0, 20.0, 25.0},

      sab { 1.19222E01, 1.14426E01, 1.00401E01, 8.05385E00, 5.90638E00, 
      3.96007E00, 2.42756E00, 1.36072E00, 6.97614E-1, 3.26730E-1, 1.40210E-1, 
      7.71234E-4, 7.07628E-4, 7.01259E-4, 7.31179E-4, 7.60155E-4, 7.90387E-4, 
      8.24931E-4, 8.59093E-4, 8.93098E-4, 1.06661E-3, 2.00558E-4, 3.38301E-7, 
      3.61904E-8, 6.3077E-11, 2.4110E-15, 5.0616E-20, 3.7100E-25, 8.39868E00, 
      8.23244E00, 7.75319E00, 6.94586E00, 5.97872E00, 4.89544E00, 3.85143E00, 
      2.88251E00, 2.07298E00, 1.41837E00, 9.32710E-1, 2.74030E-3, 1.41332E-3, 
      1.40448E-3, 1.45726E-3, 1.51636E-3, 1.57696E-3, 1.64501E-3, 1.71363E-3, 
      1.77906E-3, 2.11655E-3, 3.99875E-4, 1.36306E-6, 1.45126E-7, 5.0791E-10, 
      3.8811E-14, 1.6414E-18, 2.4920E-23, 5.89517E00, 5.85516E00, 5.68655E00, 
      5.40042E00, 5.01506E00, 4.55404E00, 4.04382E00, 3.51127E00, 2.98141E00, 
      2.47553E00, 2.01010E00, 7.55325E-2, 3.11787E-3, 2.80949E-3, 2.89853E-3, 
      3.01635E-3, 3.13811E-3, 3.27149E-3, 3.40776E-3, 3.53255E-3, 4.18102E-3, 
      7.99269E-4, 5.49535E-6, 5.86073E-7, 4.09849E-9, 6.2831E-13, 5.3903E-17, 
      1.7256E-21, 4.77865E00, 4.76124E00, 4.67504E00, 4.52376E00, 4.31386E00, 
      4.05400E00, 3.75453E00, 3.42676E00, 3.07786E00, 2.72105E00, 2.37080E00, 
      2.67105E-1, 1.08087E-2, 4.24816E-3, 4.32787E-3, 4.49995E-3, 4.68294E-3, 
      4.88013E-3, 5.08138E-3, 5.26310E-3, 6.20594E-3, 1.20111E-3, 1.24211E-5, 
      1.33362E-6, 1.39197E-8, 3.2171E-12, 4.1973E-16, 2.1043E-20, 4.10892E00, 
      4.09872E00, 4.04793E00, 3.95808E00, 3.83178E00, 3.65440E00, 3.45064E00, 
      3.22592E00, 2.98592E00, 2.72278E00, 2.45824E00, 4.84492E-1, 3.59122E-2, 
      6.23623E-3, 5.75077E-3, 5.96790E-3, 6.21127E-3, 6.47103E-3, 6.73454E-3, 
      6.97170E-3, 8.19602E-3, 1.60597E-3, 2.21489E-5, 2.39907E-6, 3.31650E-8, 
      1.0281E-11, 1.8125E-15, 1.2582E-19, 2.49247E00, 2.49319E00, 2.49389E00, 
      2.47371E00, 2.45117E00, 2.41098E00, 2.36660E00, 2.30831E00, 2.24457E00, 
      2.17099E00, 2.09128E00, 1.12574E00, 3.93892E-1, 9.67735E-2, 2.55631E-2, 
      1.55080E-2, 1.51005E-2, 1.56466E-2, 1.62433E-2, 1.67907E-2, 1.95202E-2, 
      4.10342E-3, 1.39247E-4, 1.60908E-5, 5.29253E-7, 4.2731E-10, 2.0279E-13, 
      4.1813E-17, 1.64890E00, 1.65175E00, 1.65458E00, 1.65714E00, 1.64906E00, 
      1.64093E00, 1.63243E00, 1.61377E00, 1.59532E00, 1.57644E00, 1.54819E00, 
      1.16875E00, 7.12063E-1, 3.54178E-1, 1.52300E-1, 6.59622E-2, 3.79051E-2, 
      3.14231E-2, 3.09357E-2, 3.16864E-2, 3.64635E-2, 8.50620E-3, 5.54129E-4, 
      7.14916E-5, 4.32761E-6, 7.49295E-9, 7.8951E-12, 3.9522E-15, 1.26304E00, 
      1.26596E00, 1.26898E00, 1.27216E00, 1.27281E00, 1.26878E00, 1.26485E00, 
      1.26109E00, 1.25262E00, 1.24221E00, 1.23209E00, 1.04052E00, 7.63438E-1, 
      4.90392E-1, 2.80501E-1, 1.51046E-1, 8.50637E-2, 5.72195E-2, 4.79161E-2, 
      4.59610E-2, 5.14018E-2, 1.31529E-2, 1.23250E-3, 1.75952E-4, 1.48295E-5, 
      4.11250E-8, 7.0947E-11, 6.1242E-14, 1.02854E00, 1.03134E00, 1.03423E00, 
      1.03716E00, 1.04003E00, 1.03911E00, 1.03709E00, 1.03510E00, 1.03318E00, 
      1.02863E00, 1.02211E00, 9.14765E-1, 7.37325E-1, 5.36860E-1, 3.59341E-1, 
      2.25976E-1, 1.40597E-1, 9.35593E-2, 7.11418E-2, 6.24112E-2, 6.46174E-2, 
      1.79809E-2, 2.15916E-3, 3.37909E-4, 3.55531E-5, 1.39708E-7, 3.4654E-10, 
      4.4497E-13, 4.67971E-1, 4.69718E-1, 4.71191E-1, 4.72620E-1, 4.74307E-1, 
      4.76151E-1, 4.77962E-1, 4.79340E-1, 4.79879E-1, 4.80325E-1, 4.80818E-1, 
      4.78303E-1, 4.56328E-1, 4.22714E-1, 3.76759E-1, 3.28624E-1, 2.79851E-1, 
      2.36785E-1, 1.99551E-1, 1.71042E-1, 1.16902E-1, 4.80050E-2, 1.19814E-2, 
      2.83678E-3, 5.67767E-4, 7.45367E-6, 6.34067E-8, 3.0367E-10, 2.12962E-1, 
      2.14002E-1, 2.15155E-1, 2.16043E-1, 2.16794E-1, 2.17565E-1, 2.18358E-1, 
      2.19172E-1, 2.19834E-1, 2.20573E-1, 2.21338E-1, 2.28527E-1, 2.30922E-1, 
      2.31912E-1, 2.28484E-1, 2.23803E-1, 2.16477E-1, 2.08341E-1, 1.99454E-1, 
      1.90369E-1, 1.45786E-1, 8.88424E-2, 3.71167E-2, 1.34313E-2, 4.24235E-3, 
      1.54442E-4, 3.61071E-6, 4.94098E-8, 1.22316E-1, 1.22894E-1, 1.23558E-1, 
      1.24324E-1, 1.25106E-1, 1.25701E-1, 1.26179E-1, 1.26682E-1, 1.27067E-1, 
      1.27468E-1, 1.27954E-1, 1.33479E-1, 1.37814E-1, 1.41464E-1, 1.44185E-1, 
      1.45690E-1, 1.46945E-1, 1.47364E-1, 1.47666E-1, 1.47167E-1, 1.39640E-1, 
      1.09656E-1, 6.21654E-2, 2.97587E-2, 1.23539E-2, 8.54115E-4, 3.73825E-5, 
      9.70721E-7, 7.86262E-2, 7.90164E-2, 7.94156E-2, 7.98240E-2, 8.03053E-2, 
      8.08288E-2, 8.12942E-2, 8.16809E-2, 8.20454E-2, 8.23623E-2, 8.26832E-2, 
      8.65317E-2, 8.97016E-2, 9.31749E-2, 9.60392E-2, 9.90652E-2, 1.01521E-1, 
      1.03591E-1, 1.05794E-1, 1.07824E-1, 1.18720E-1, 1.13495E-1, 8.03171E-2, 
      4.73309E-2, 2.40126E-2, 2.67133E-3, 1.85773E-4, 7.70570E-6, 5.38110E-2, 
      5.40757E-2, 5.43447E-2, 5.46181E-2, 5.48958E-2, 5.51834E-2, 5.55871E-2, 
      5.59938E-2, 5.62991E-2, 5.64973E-2, 5.67019E-2, 5.93824E-2, 6.17773E-2, 
      6.45718E-2, 6.71294E-2, 6.93127E-2, 7.19010E-2, 7.42549E-2, 7.61853E-2, 
      7.85081E-2, 9.55268E-2, 1.06826E-1, 8.98468E-2, 6.24135E-2, 3.71358E-2, 
      6.04519E-3, 6.08705E-4, 3.66195E-5, 2.38647E-2, 2.39834E-2, 2.41057E-2, 
      2.42390E-2, 2.43734E-2, 2.45090E-2, 2.46451E-2, 2.47822E-2, 2.49110E-2, 
      2.50441E-2, 2.51807E-2, 2.63939E-2, 2.76416E-2, 2.88719E-2, 3.02718E-2, 
      3.16042E-2, 3.28423E-2, 3.42612E-2, 3.56453E-2, 3.70513E-2, 5.10213E-2, 
      7.48974E-2, 8.51303E-2, 7.97293E-2, 6.36667E-2, 2.11653E-2, 4.27511E-3, 
      5.15329E-4, 1.16097E-2, 1.16677E-2, 1.17284E-2, 1.17897E-2, 1.18513E-2, 
      1.19132E-2, 1.19753E-2, 1.20374E-2, 1.20982E-2, 1.21581E-2, 1.22185E-2, 
      1.28305E-2, 1.34632E-2, 1.40865E-2, 1.47695E-2, 1.54604E-2, 1.61699E-2, 
      1.68653E-2, 1.76069E-2, 1.83788E-2, 2.67444E-2, 4.60227E-2, 6.33204E-2, 
      7.26100E-2, 7.12479E-2, 3.96974E-2, 1.33332E-2, 2.66750E-3, 5.89844E-3, 
      5.92797E-3, 5.95762E-3, 5.98741E-3, 6.01753E-3, 6.04805E-3, 6.07869E-3, 
      6.10943E-3, 6.14005E-3, 6.17065E-3, 6.20139E-3, 6.51531E-3, 6.83930E-3, 
      7.16968E-3, 7.51302E-3, 7.87430E-3, 8.24531E-3, 8.62657E-3, 9.01387E-3, 
      9.41929E-3, 1.41296E-2, 2.68360E-2, 4.19896E-2, 5.56388E-2, 6.36015E-2, 
      5.27244E-2, 2.63599E-2, 7.85109E-3, 3.07259E-3, 3.08799E-3, 3.10347E-3, 
      3.11902E-3, 3.13461E-3, 3.15026E-3, 3.16600E-3, 3.18185E-3, 3.19770E-3, 
      3.21360E-3, 3.22957E-3, 3.39368E-3, 3.56298E-3, 3.73781E-3, 3.91997E-3, 
      4.10960E-3, 4.30739E-3, 4.51051E-3, 4.71992E-3, 4.93732E-3, 7.54306E-3, 
      1.52908E-2, 2.61961E-2, 3.86255E-2, 4.96294E-2, 5.64350E-2, 3.89692E-2, 
      1.60735E-2, 1.62762E-3, 1.63578E-3, 1.64399E-3, 1.65223E-3, 1.66051E-3, 
      1.66881E-3, 1.67711E-3, 1.68544E-3, 1.69381E-3, 1.70223E-3, 1.71069E-3, 
      1.79722E-3, 1.88729E-3, 1.98114E-3, 2.07871E-3, 2.18029E-3, 2.28553E-3, 
      2.39490E-3, 2.50865E-3, 2.62646E-3, 4.06319E-3, 8.62231E-3, 1.57874E-2, 
      2.52482E-2, 3.55534E-2, 5.22565E-2, 4.71874E-2, 2.55802E-2, 8.72309E-4, 
      8.76681E-4, 8.81072E-4, 8.85485E-4, 8.89920E-4, 8.94377E-4, 8.98856E-4, 
      9.03357E-4, 9.07882E-4, 9.12424E-4, 9.16967E-4, 9.63346E-4, 1.01168E-3, 
      1.06231E-3, 1.11506E-3, 1.16972E-3, 1.22688E-3, 1.28611E-3, 1.34767E-3, 
      1.41202E-3, 2.20405E-3, 4.83705E-3, 9.31501E-3, 1.58711E-2, 2.40493E-2, 
      4.36911E-2, 4.94837E-2, 3.38751E-2, 4.71648E-4, 4.74014E-4, 4.76390E-4, 
      4.78778E-4, 4.81177E-4, 4.83588E-4, 4.86003E-4, 4.88429E-4, 4.90868E-4, 
      4.93318E-4, 4.95780E-4, 5.20952E-4, 5.47189E-4, 5.74643E-4, 6.03219E-4, 
      6.32998E-4, 6.64071E-4, 6.96423E-4, 7.30133E-4, 7.65281E-4, 1.20263E-3, 
      2.70794E-3, 5.42421E-3, 9.72025E-3, 1.56345E-2, 3.39057E-2, 4.66052E-2, 
      3.90322E-2, 2.56817E-4, 2.58096E-4, 2.59382E-4, 2.60674E-4, 2.61973E-4, 
      2.63277E-4, 2.64587E-4, 2.65904E-4, 2.67227E-4, 2.68557E-4, 2.69894E-4, 
      2.83617E-4, 2.98021E-4, 3.12951E-4, 3.28518E-4, 3.44835E-4, 3.61868E-4, 
      3.79562E-4, 3.98097E-4, 4.17510E-4, 6.59481E-4, 1.51539E-3, 3.13253E-3, 
      5.84823E-3, 9.88108E-3, 2.48891E-2, 4.04424E-2, 4.03814E-2},

      correct { 1.19160E01, 1.14366E01, 1.00349E01, 8.04966E00, 5.90330E00, 
      3.95801E00, 2.42629E00, 1.36001E00, 6.97250E-1, 3.26559E-1, 1.40137E-1, 
      7.70832E-4, 7.07259E-4, 7.00894E-4, 7.30798E-4, 7.59759E-4, 7.89976E-4, 
      8.24502E-4, 8.58646E-4, 8.92633E-4, 1.06605E-3, 2.00496E-4, 9.75623E-4, 
      3.22805E-7, 4.86804E-8, 2.16534E-7, 2.35699E-9, 3.7054E-13, 8.38993E00, 
      8.22387E00, 7.74511E00, 6.93862E00, 5.97250E00, 4.89034E00, 3.84742E00, 
      2.87951E00, 2.07082E00, 1.41689E00, 9.31738E-1, 2.73744E-3, 1.41185E-3, 
      1.40302E-3, 1.45574E-3, 1.51478E-3, 1.57532E-3, 1.64330E-3, 1.71184E-3, 
      1.77720E-3, 2.11435E-3, 3.99627E-4, 2.57822E-3, 1.28210E-6, 1.97846E-7, 
      8.62052E-7, 5.17659E-8, 4.4025E-12, 5.88289E00, 5.84296E00, 5.67471E00, 
      5.38918E00, 5.00462E00, 4.54456E00, 4.03540E00, 3.50396E00, 2.97520E00, 
      2.47038E00, 2.00592E00, 7.53752E-2, 3.11138E-3, 2.80364E-3, 2.89250E-3, 
      3.01007E-3, 3.13158E-3, 3.26468E-3, 3.40067E-3, 3.52520E-3, 4.17233E-3, 
      7.98268E-4, 4.95285E-3, 5.07412E-6, 8.06729E-7, 3.41911E-6, 8.87400E-7, 
      5.2435E-11, 4.76373E00, 4.74638E00, 4.66044E00, 4.50964E00, 4.30039E00, 
      4.04135E00, 3.74281E00, 3.41606E00, 3.06825E00, 2.71256E00, 2.36340E00, 
      2.66271E-1, 1.07749E-2, 4.23490E-3, 4.31436E-3, 4.48590E-3, 4.66832E-3, 
      4.86490E-3, 5.06552E-3, 5.24667E-3, 6.18661E-3, 1.19883E-3, 6.68663E-3, 
      1.13175E-5, 1.84080E-6, 7.63231E-6, 3.17232E-6, 2.2370E-10, 4.09183E00, 
      4.08167E00, 4.03109E00, 3.94161E00, 3.81584E00, 3.63920E00, 3.43629E00, 
      3.21250E00, 2.97350E00, 2.71145E00, 2.44801E00, 4.82476E-1, 3.57628E-2, 
      6.21029E-3, 5.72685E-3, 5.94308E-3, 6.18544E-3, 6.44412E-3, 6.70653E-3, 
      6.94270E-3, 8.16200E-3, 1.60188E-3, 8.07727E-3, 1.99645E-5, 3.31157E-6, 
      1.34653E-5, 6.80723E-6, 1.05905E-9, 2.46663E00, 2.46734E00, 2.46804E00, 
      2.44806E00, 2.42575E00, 2.38598E00, 2.34207E00, 2.28438E00, 2.22130E00, 
      2.14848E00, 2.06960E00, 1.11406E00, 3.89809E-1, 9.57702E-2, 2.52981E-2, 
      1.53473E-2, 1.49440E-2, 1.54844E-2, 1.60749E-2, 1.66167E-2, 1.93183E-2, 
      4.07624E-3, 1.34791E-2, 1.20114E-4, 2.18497E-5, 8.05937E-5, 4.76384E-5, 
      2.38789E-8, 1.61489E00, 1.61768E00, 1.62045E00, 1.62295E00, 1.61505E00, 
      1.60708E00, 1.59876E00, 1.58048E00, 1.56241E00, 1.54392E00, 1.51626E00, 
      1.14465E00, 6.97375E-1, 3.46872E-1, 1.49158E-1, 6.46017E-2, 3.71233E-2, 
      3.07750E-2, 3.02977E-2, 3.10330E-2, 3.57132E-2, 8.38758E-3, 1.85794E-2, 
      4.56304E-4, 9.33638E-5, 3.01459E-4, 1.55947E-4, 2.61768E-7, 1.22417E00, 
      1.22700E00, 1.22992E00, 1.23300E00, 1.23363E00, 1.22973E00, 1.22592E00, 
      1.22227E00, 1.21407E00, 1.20397E00, 1.19416E00, 1.00849E00, 7.39940E-1, 
      4.75298E-1, 2.71868E-1, 1.46397E-1, 8.24457E-2, 5.54587E-2, 4.64416E-2, 
      4.45468E-2, 4.98240E-2, 1.28671E-2, 2.19086E-2, 9.80885E-4, 2.21144E-4, 
      6.49308E-4, 2.92135E-4, 1.08572E-6, 9.86546E-1, 9.89240E-1, 9.92008E-1, 
      9.94816E-1, 9.97573E-1, 9.96687E-1, 9.94755E-1, 9.92847E-1, 9.91001E-1, 
      9.86633E-1, 9.80384E-1, 8.77417E-1, 7.07222E-1, 5.14941E-1, 3.44670E-1, 
      2.16750E-1, 1.34857E-1, 8.97402E-2, 6.82381E-2, 5.98641E-2, 6.19870E-2, 
      1.74452E-2, 2.44627E-2, 1.67243E-3, 4.09972E-4, 1.15943E-3, 4.46623E-4, 
      3.01028E-6, 4.21664E-1, 4.23237E-1, 4.24565E-1, 4.25852E-1, 4.27372E-1, 
      4.29034E-1, 4.30665E-1, 4.31908E-1, 4.32393E-1, 4.32794E-1, 4.33239E-1, 
      4.30974E-1, 4.11173E-1, 3.80886E-1, 3.39480E-1, 2.96108E-1, 2.52163E-1, 
      2.13360E-1, 1.79811E-1, 1.54125E-1, 1.05382E-1, 4.41159E-2, 3.49296E-2, 
      8.34395E-3, 2.92398E-3, 7.34235E-3, 1.66118E-3, 7.47140E-5, 1.72918E-1, 
      1.73763E-1, 1.74699E-1, 1.75420E-1, 1.76030E-1, 1.76656E-1, 1.77300E-1, 
      1.77961E-1, 1.78499E-1, 1.79098E-1, 1.79720E-1, 1.85558E-1, 1.87504E-1, 
      1.88309E-1, 1.85528E-1, 1.81730E-1, 1.75785E-1, 1.69183E-1, 1.61972E-1, 
      1.54601E-1, 1.18523E-1, 7.40554E-2, 5.00074E-2, 2.40571E-2, 1.16243E-2, 
      1.58750E-2, 4.84991E-3, 7.15007E-4, 8.95226E-2, 8.99459E-2, 9.04316E-2, 
      9.09920E-2, 9.15645E-2, 9.19998E-2, 9.23503E-2, 9.27184E-2, 9.30003E-2, 
      9.32934E-2, 9.36496E-2, 9.76937E-2, 1.00867E-1, 1.03541E-1, 1.05536E-1, 
      1.06642E-1, 1.07565E-1, 1.07878E-1, 1.08106E-1, 1.07749E-1, 1.02428E-1, 
      8.26883E-2, 6.10011E-2, 3.89297E-2, 2.28159E-2, 1.96255E-2, 9.22553E-3, 
      2.35463E-3, 5.18930E-2, 5.21506E-2, 5.24140E-2, 5.26836E-2, 5.30012E-2, 
      5.33466E-2, 5.36537E-2, 5.39090E-2, 5.41496E-2, 5.43588E-2, 5.45706E-2, 
      5.71111E-2, 5.92049E-2, 6.14993E-2, 6.33929E-2, 6.53939E-2, 6.70203E-2, 
      6.83929E-2, 6.98545E-2, 7.12031E-2, 7.86043E-2, 7.73225E-2, 6.49505E-2, 
      4.88433E-2, 3.30794E-2, 2.23391E-2, 1.39487E-2, 5.00363E-3, 3.20426E-2, 
      3.22002E-2, 3.23604E-2, 3.25232E-2, 3.26886E-2, 3.28598E-2, 3.31000E-2, 
      3.33419E-2, 3.35237E-2, 3.36419E-2, 3.37638E-2, 3.53604E-2, 3.67878E-2, 
      3.84534E-2, 3.99789E-2, 4.12829E-2, 4.28284E-2, 4.42357E-2, 4.53925E-2, 
      4.67837E-2, 5.71139E-2, 6.57232E-2, 6.26148E-2, 5.31570E-2, 4.04452E-2, 
      2.52766E-2, 1.83206E-2, 8.36286E-3, 1.10127E-2, 1.10674E-2, 1.11239E-2, 
      1.11853E-2, 1.12473E-2, 1.13099E-2, 1.13726E-2, 1.14358E-2, 1.14952E-2, 
      1.15565E-2, 1.16195E-2, 1.21795E-2, 1.27557E-2, 1.33244E-2, 1.39714E-2, 
      1.45880E-2, 1.51622E-2, 1.58198E-2, 1.64623E-2, 1.71157E-2, 2.36782E-2, 
      3.57089E-2, 4.37513E-2, 4.65658E-2, 4.42092E-2, 3.22369E-2, 2.62335E-2, 
      1.73174E-2, 4.16314E-3, 4.18395E-3, 4.20571E-3, 4.22768E-3, 4.24976E-3, 
      4.27196E-3, 4.29416E-3, 4.31641E-3, 4.33818E-3, 4.35966E-3, 4.38128E-3, 
      4.60064E-3, 4.82758E-3, 5.05152E-3, 5.29684E-3, 5.54528E-3, 5.80068E-3, 
      6.05153E-3, 6.31916E-3, 6.59792E-3, 9.65182E-3, 1.70262E-2, 2.47608E-2, 
      3.11145E-2, 3.47414E-2, 3.35524E-2, 2.97839E-2, 2.38315E-2, 1.64672E-3, 
      1.65497E-3, 1.66324E-3, 1.67156E-3, 1.67997E-3, 1.68849E-3, 1.69703E-3, 
      1.70560E-3, 1.71413E-3, 1.72266E-3, 1.73123E-3, 1.81877E-3, 1.90921E-3, 
      2.00154E-3, 2.09758E-3, 2.19869E-3, 2.30266E-3, 2.40963E-3, 2.51848E-3, 
      2.63252E-3, 3.97057E-3, 7.71400E-3, 1.26253E-2, 1.79700E-2, 2.27455E-2, 
      2.87564E-2, 2.92576E-2, 2.65997E-2, 6.68664E-4, 6.72016E-4, 6.75385E-4, 
      6.78768E-4, 6.82162E-4, 6.85569E-4, 6.88990E-4, 6.92431E-4, 6.95874E-4, 
      6.99326E-4, 7.02796E-4, 7.38450E-4, 7.75272E-4, 8.13339E-4, 8.53035E-4, 
      8.94395E-4, 9.37573E-4, 9.81983E-4, 1.02783E-3, 1.07548E-3, 1.65200E-3, 
      3.41896E-3, 6.08638E-3, 9.52266E-3, 1.33078E-2, 2.12863E-2, 2.54361E-2, 
      2.59464E-2, 2.76324E-4, 2.77710E-4, 2.79103E-4, 2.80502E-4, 2.81908E-4, 
      2.83318E-4, 2.84725E-4, 2.86136E-4, 2.87554E-4, 2.88980E-4, 2.90412E-4, 
      3.05075E-4, 3.20353E-4, 3.36282E-4, 3.52862E-4, 3.70138E-4, 3.88059E-4, 
      4.06704E-4, 4.26118E-4, 4.46252E-4, 6.94005E-4, 1.50097E-3, 2.84224E-3, 
      4.78392E-3, 7.23702E-3, 1.41546E-2, 1.99039E-2, 2.28200E-2, 1.15589E-4, 
      1.16169E-4, 1.16751E-4, 1.17336E-4, 1.17923E-4, 1.18514E-4, 1.19106E-4, 
      1.19701E-4, 1.20299E-4, 1.20899E-4, 1.21499E-4, 1.27632E-4, 1.34029E-4, 
      1.40734E-4, 1.47726E-4, 1.54981E-4, 1.62574E-4, 1.70452E-4, 1.78650E-4, 
      1.87227E-4, 2.93735E-4, 6.56019E-4, 1.30180E-3, 2.32068E-3, 3.74404E-3, 
      8.69396E-3, 1.42586E-2, 1.83987E-2, 4.87972E-5, 4.90420E-5, 4.92879E-5, 
      4.95351E-5, 4.97834E-5, 5.00329E-5, 5.02822E-5, 5.05324E-5, 5.07840E-5, 
      5.10367E-5, 5.12906E-5, 5.38884E-5, 5.65988E-5, 5.94366E-5, 6.23936E-5, 
      6.54781E-5, 6.86998E-5, 7.20580E-5, 7.55610E-5, 7.92175E-5, 1.25100E-4, 
      2.86286E-4, 5.89197E-4, 1.09978E-3, 1.87067E-3, 5.03033E-3, 9.49887E-3, 
      1.37652E-2, 2.07505E-5, 2.08539E-5, 2.09579E-5, 2.10624E-5, 2.11674E-5, 
      2.12729E-5, 2.13786E-5, 2.14846E-5, 2.15912E-5, 2.16984E-5, 2.18061E-5, 
      2.29119E-5, 2.40731E-5, 2.52784E-5, 2.65364E-5, 2.78558E-5, 2.92344E-5, 
      3.06685E-5, 3.21720E-5, 3.37478E-5, 5.35591E-5, 1.24939E-4, 2.64715E-4, 
      5.12970E-4, 9.11703E-4, 2.78109E-3, 5.96512E-3, 9.66009E-3 };

      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, temp, 
                  oscEnergiesWeights, effectiveTemp, sab );
      THEN( "scattering law matrix is correctly changed" ){
         REQUIRE(ranges::equal(sab,correct,equal));
      } // THEN
    } // WHEN
    WHEN( "alpha and beta values are scaled (lat=1)" ){
      std::vector<double>
      alpha {0.01, 0.08, 0.40, 0.80, 4.0, 6.0, 20., 30, 40.0, 50},
      beta  {0.00, 0.04, 0.06, 0.08, 0.1, 0.3, 0.5, 0.9, 2.0, 6.0, 10, 20},
      sab = { 1.39300E+1, 8.31737E+0, 4.33327E+0, 1.72907E+0, 5.31649E-1, 
        5.44280E-4, 4.99581E-4, 5.73572E-4, 7.32629E-4, 1.18469E-4, 4.33780E-8, 
        5.8434E-16, 4.83282E+0, 4.57996E+0, 4.25367E+0, 3.81799E+0, 3.31165E+0, 
        1.35701E-1, 4.22317E-3, 4.50697E-3, 5.67766E-3, 9.42199E-4, 2.85498E-6, 
        2.6454E-12, 1.99384E+0, 1.99446E+0, 1.97516E+0, 1.94448E+0, 1.90383E+0, 
        1.06653E+0, 3.20660E-1, 2.45166E-2, 2.58739E-2, 4.87027E-3, 7.75783E-5, 
        2.00293E-9, 1.28258E+0, 1.29134E+0, 1.29017E+0, 1.28422E+0, 1.27762E+0, 
        9.97146E-1, 5.71174E-1, 1.00587E-1, 4.69977E-2, 1.01645E-2, 3.34362E-4, 
        3.71882E-8, 3.10519E-1, 3.13549E-1, 3.15279E-1, 3.17153E-1, 3.19297E-1, 
        3.26232E-1, 3.18139E-1, 2.68403E-1, 1.40426E-1, 5.64745E-2, 1.02841E-2, 
        3.99469E-5, 1.88287E-1, 1.90836E-1, 1.91814E-1, 1.92833E-1, 1.93900E-1, 
        2.03716E-1, 2.07823E-1, 2.01375E-1, 1.51823E-1, 7.96267E-2, 2.24716E-2, 
        2.28227E-4, 2.48679E-2, 2.52349E-2, 2.54277E-2, 2.56389E-2, 2.58515E-2, 
        2.77421E-2, 2.95846E-2, 3.35846E-2, 4.46011E-2, 7.76412E-2, 8.08930E-2, 
        1.76304E-2, 8.71638E-3, 8.84808E-3, 8.91473E-3, 8.98192E-3, 9.04964E-3, 
        9.72164E-3, 1.04265E-2, 1.19271E-2, 1.66844E-2, 3.97662E-2, 6.12275E-2, 
        3.83700E-2, 3.30664E-3, 3.35554E-3, 3.38034E-3, 3.40538E-3, 3.43050E-3, 
        3.68809E-3, 3.95922E-3, 4.54709E-3, 6.49196E-3, 1.81904E-2, 3.53655E-2, 
        4.51282E-2, 1.30443E-3, 1.32372E-3, 1.33345E-3, 1.34317E-3, 1.35297E-3, 
        1.45449E-3, 1.56242E-3, 1.79816E-3, 2.59674E-3, 7.99775E-3, 1.80875E-2, 
        3.83645E-2},
      correct  { 1.39247E+1, 8.31418E+0, 4.33161E+0, 1.72841E+0, 5.31445E-1, 
        5.44072E-4, 4.99390E-4, 5.73353E-4, 7.32349E-4, 1.18455E-4, 1.90979E-7, 
        1.03644E-7, 4.81805E+0, 4.56597E+0, 4.24067E+0, 3.80633E+0, 3.30154E+0, 
        1.35286E-1, 4.21027E-3, 4.49320E-3, 5.66033E-3, 9.41226E-4, 1.19856E-5, 
        6.49077E-6, 1.96357E+0, 1.96418E+0, 1.94517E+0, 1.91496E+0, 1.87493E+0, 
        1.05034E+0, 3.15792E-1, 2.41444E-2, 2.54815E-2, 4.83935E-3, 2.85548E-4, 
        1.71243E-4, 1.24393E+0, 1.25243E+0, 1.25129E+0, 1.24552E+0, 1.23912E+0, 
        9.67099E-1, 5.53963E-1, 9.75569E-2, 4.55832E-2, 1.00126E-2, 1.13174E-3, 
        1.25942E-3, 2.66483E-1, 2.69083E-1, 2.70568E-1, 2.72176E-1, 2.74016E-1, 
        2.79968E-1, 2.73025E-1, 2.30347E-1, 1.20553E-1, 5.05390E-2, 1.93912E-2, 
        1.55770E-2, 1.49705E-1, 1.51732E-1, 1.52509E-1, 1.53320E-1, 1.54169E-1, 
        1.61974E-1, 1.65242E-1, 1.60124E-1, 1.20787E-1, 6.64410E-2, 3.32096E-2, 
        1.81356E-2, 1.16575E-2, 1.18295E-2, 1.19199E-2, 1.20188E-2, 1.21184E-2, 
        1.30054E-2, 1.38710E-2, 1.57526E-2, 2.09744E-2, 3.80279E-2, 4.60660E-2, 
        2.78409E-2, 2.82117E-3, 2.86378E-3, 2.88535E-3, 2.90709E-3, 2.92901E-3, 
        3.14676E-3, 3.37549E-3, 3.86351E-3, 5.42303E-3, 1.33821E-2, 2.26963E-2, 
        2.68342E-2, 7.41522E-4, 7.52491E-4, 7.58053E-4, 7.63668E-4, 7.69304E-4, 
        8.27141E-4, 8.88117E-4, 1.02065E-3, 1.46265E-3, 4.22139E-3, 8.83731E-3, 
        1.79316E-2, 2.03062E-4, 2.06066E-4, 2.07581E-4, 2.09096E-4, 2.10622E-4, 
        2.26451E-4, 2.43303E-4, 2.80205E-4, 4.06179E-4, 1.28286E-3, 3.08555E-3, 
        9.46688E-3};
      double sc = 0.7339860234, scaling = 0.7339860234;
      for (auto& a : alpha){ a*= scaling; }
      for (auto& b : beta ){ b*= sc;      }


      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, temp, 
              oscEnergiesWeights, effectiveTemp, sab );

      THEN( "scattering law matrix is correctly changed" ){
         REQUIRE(ranges::equal(sab,correct,equal));
      } // THEN
    } // WHEN
  } // GIVEN
} // TEST CASE 























TEST_CASE( "Discrete oscillator treatment (old tests)" ){
  double temp  = 200.0012;
  GIVEN( "two oscillators" ){
    WHEN( "alpha and beta values are slightly small" ){
      std::vector<double> osc_energies{0.035, 0.05}, osc_weights{0.2, 0.8},
        alpha{0.1, 0.2, 0.4, 0.8, 1.6}, beta{0.10, 0.15, 0.30, 0.60, 1.20};
      double t_eff = 81178.935219;

      std::vector<double> sym_sab (alpha.size()*beta.size(),0.0);
      for ( size_t i = 0; i < sym_sab.size(); ++i ){ sym_sab[i] = i+1; }

      double lambda_s = 2.2941534E-3, tbeta = 2, twt = 0.3;

    auto oscEnergiesWeights = ranges::view::zip(osc_energies,osc_weights);

      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, 
              temp, oscEnergiesWeights, t_eff, sym_sab );

      std::vector<double> correctSymSab {0.9575582, 1.914953, 2.872356, 
        3.829659, 4.807700, 5.501617, 6.418366, 7.335187, 8.251386, 9.253313, 
        9.256188, 10.09717, 10.93859, 11.77685, 12.85839, 11.36789, 12.07716, 
        12.78842, 13.48618, 14.74575, 10.74546, 11.25461, 11.77105, 12.24116, 
        13.75305};

      THEN( "scattering law matrix is correctly changed" ){
        REQUIRE(ranges::equal(sym_sab,correctSymSab,equal));
      } // THEN
    } // WHEN

    WHEN( "alpha and beta values are slightly larger" ){
      std::vector<double> osc_energies{0.035, 0.05}, osc_weights{0.2, 0.8},
        alpha{2.1, 4.2, 6.4, 8.8, 10.6}, beta{1.10, 2.15, 3.30, 4.60, 5.20};
      double t_eff = 81178.935219;

      std::vector<double> sym_sab (alpha.size()*beta.size(),0.0);
      for ( size_t i = 0; i < sym_sab.size(); ++i ){ sym_sab[i] = i+1; }


      double lambda_s = 2.2941534E-3, tbeta = 2.0, twt = 0.3;

    auto oscEnergiesWeights = ranges::view::zip(osc_energies,osc_weights);

      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, 
              temp, oscEnergiesWeights, t_eff, sym_sab );

      std::vector<double> correctSymSab {0.7125247, 1.321226, 1.681776, 2.414681,
        3.053423, 2.268664, 3.482821, 3.839385, 5.171756, 5.720471, 2.353392, 
        3.703027, 4.256384, 6.170473, 6.777184, 1.859812, 2.977392, 3.629533, 
        5.625912, 6.256402, 1.551690, 2.499096, 3.173748, 5.098916, 5.735165};

      THEN( "scattering law matrix is correctly changed" ){
        REQUIRE(ranges::equal(sym_sab,correctSymSab,equal));
      } // THEN
    } // WHEN

    WHEN( "oscillator weights add up to 1" ){
      std::vector<double> osc_energies{0.1, 0.5}, osc_weights{0.4, 0.6},
        alpha{2.1, 4.2, 6.4, 8.8, 10.6}, beta{1.10, 2.15, 3.30, 4.60, 5.20};
      double t_eff = 81178.935219;

      std::vector<double> sym_sab (alpha.size()*beta.size(),0.0);
      for ( size_t i = 0; i < sym_sab.size(); ++i ){ sym_sab[i] = i+1; }

      double lambda_s = 2.2941534E-3, tbeta = 2.0, twt = 0.3;

    auto oscEnergiesWeights = ranges::view::zip(osc_energies,osc_weights);

      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, 
              temp, oscEnergiesWeights, t_eff, sym_sab );

      std::vector<double> correctSymSab {0.8323044, 1.665653, 2.505608, 
        3.415392, 4.259058, 4.128312, 4.839682, 5.601205, 7.150784, 8.046933, 
        6.217680, 6.833610, 7.560954, 10.11139, 11.17169, 7.300605, 7.836529, 
        8.549848, 12.21146, 13.47660, 8.160909, 8.655053, 9.386204, 14.08904, 
        15.57309};

      THEN( "scattering law matrix is correctly changed" ){
        REQUIRE(ranges::equal(sym_sab,correctSymSab,equal));
      } // THEN
    } // WHEN
    WHEN( "oscillator weights don't add up to 1" ){
      std::vector<double> osc_energies{0.1, 0.5}, osc_weights{0.4, 0.5},
        alpha{2.1, 4.2, 6.4, 8.8, 10.6}, beta{1.10, 2.15, 3.30, 4.60, 5.20};
      double t_eff = 81178.935219;

      std::vector<double> sym_sab (alpha.size()*beta.size(),0.0);
      for ( size_t i = 0; i < sym_sab.size(); ++i ){ sym_sab[i] = i+1; }

      double lambda_s = 2.2941534E-3, tbeta = 2.0, twt = 0.3;

    auto oscEnergiesWeights = ranges::view::zip(osc_energies,osc_weights);
      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, 
              temp, oscEnergiesWeights, t_eff, sym_sab );

      std::vector<double> correctSymSab {0.83835105, 1.6777545, 2.5238118, 
        3.4402043, 4.2899999, 4.1885133, 4.9102577, 5.6828846, 7.2550608, 
        8.1642777, 6.3563702, 6.9860384, 7.7296068, 10.336937, 11.420888, 
        7.5254485, 8.0778780, 8.8131661, 12.587555, 13.891658, 8.4646043, 
        8.9771372, 9.7354967, 14.613342, 16.15261};

      THEN( "scattering law matrix is correctly changed" ){
        REQUIRE(ranges::equal(sym_sab,correctSymSab,equal));
      } // THEN
    } // WHEN
  } // GIVEN

  GIVEN( "3 or more oscillators" ){
    WHEN( "3 oscillator weights add up to 1" ){
      std::vector<double> osc_energies{0.1, 0.2, 0.3}, osc_weights{0.2, 0.3, 0.5},
        alpha{2.1, 4.2, 6.4, 8.8, 10.6}, beta{1.10, 2.15, 3.30, 4.60, 5.20};
      double t_eff = 81178.935219;

      std::vector<double> sym_sab (alpha.size()*beta.size(),0.0);
      for ( size_t i = 0; i < sym_sab.size(); ++i ){ sym_sab[i] = i+1; }

      double lambda_s = 2.2941534E-3, tbeta = 2.0, twt = 0.3;

    auto oscEnergiesWeights = ranges::view::zip(osc_energies,osc_weights);
      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, 
              temp, oscEnergiesWeights, t_eff, sym_sab );

      std::vector<double> correctSymSab {0.8313648, 1.663253, 2.498450, 
        3.368623, 4.205679, 4.132822, 4.833326, 5.558994, 6.680081, 7.473300, 
        6.230090, 6.821939, 7.469791, 9.034025, 9.849238, 7.320823, 7.818479, 
        8.405469, 10.47719, 11.34208, 8.187767, 8.630877, 9.193504, 11.75896, 
        12.70126}; 

      THEN( "scattering law matrix is correctly changed" ){
        REQUIRE(ranges::equal(sym_sab,correctSymSab,equal));
      } // THEN
    } // WHEN
    
    WHEN( "5 oscillators where weights don't add up to 1" ){
      std::vector<double> osc_energies{1,2,3,4,5}, osc_weights{0.5,0.4,0.3,0.2,0.1},
        alpha{2.1, 4.2, 6.4, 8.8, 10.6}, beta{1.10, 2.15, 3.30, 4.60, 5.20};
      double t_eff = 81178.935219;

      std::vector<double> sym_sab (alpha.size()*beta.size(),0.0);
      for ( size_t i = 0; i < sym_sab.size(); ++i ){ sym_sab[i] = i+1; }

      double lambda_s = 2.2941534E-3, tbeta = 2.0, twt = 0.3;

    auto oscEnergiesWeights = ranges::view::zip(osc_energies,osc_weights);
      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, 
              temp, oscEnergiesWeights, t_eff, sym_sab );

      std::vector<double> correctSymSab {0.9690026, 1.938005, 2.907007, 
        3.876010, 4.845013, 5.633796, 6.572763, 7.511729, 8.450695, 9.389661, 
        9.993471, 10.90196, 11.81046, 12.71896, 13.62746, 14.02216, 14.89855, 
        15.77493, 16.65132, 17.52770, 17.91401, 18.76706, 19.62011, 20.47316, 
        21.32620}; 

      THEN( "scattering law matrix is correctly changed" ){
        REQUIRE(ranges::equal(sym_sab,correctSymSab,equal));
      } // THEN
    } // WHEN
  } // GIVEN


  GIVEN( "dummy test" ){
    double lambda_s = 0.26460498561058793, temp = 296.0,
    twt = 0.0, tbeta = 0.5, effectiveTemp = 572.61028482016525;
    std::vector<double> 
      alpha = {0.1, 1.0, 10},
      beta  = {0.0, 5.0, 50.0},
      oscEnergies { 0.02, 0.05 },
      oscWeights  { 0.2, 0.3 },
      sab { 6.280552E-03, 1.484229E-03, 1.781044E-29, 5.104245E-02,
            1.522726E-02, 4.744955E-20, 7.013400E-02, 8.693012E-02,
            8.708265E-10 },
      sabCorrect { 1.434008E-02, 1.499370E-03, 0.000000E+00, 6.461972E-02,
                   1.711452E-02, 1.054588E-19, 9.475565E-03, 4.652109E-02,
                   2.206891E-08 };

    auto oscEnergiesWeights = ranges::view::zip(oscEnergies,oscWeights);

    discreteOscillators( lambda_s, twt, tbeta, alpha, beta, temp, oscEnergiesWeights, effectiveTemp, sab );
    THEN( "scattering law matrix is correctly changed" ){
      REQUIRE(ranges::equal(sab,sabCorrect,equal));
    } // THEN
  } // GIVEN



  GIVEN( "simple water example" ){
    double lambda_s = 0.23520419644942447, temp = 296.0,
    twt = 0.055556, tbeta = 0.444444, effectiveTemp = 541.87556285322705;
    std::vector<double> oscEnergies { 0.205, 0.48 }, 
                        oscWeights  { 0.166667, 0.333333 };
    auto oscEnergiesWeights = ranges::view::zip(oscEnergies,oscWeights);

    WHEN( "alpha and beta valeus are moderately sized" ){
      std::vector<double> 
        alpha = {0.01, 0.04, 0.08, 0.10, 0.40, 0.80, 1.00, 4.00, 8.00, 10.00},
        beta  = {0.00, 1.00, 2.00, 3.00, 4.00},
        sab { 1.193900E+01, 3.646870E-4, 4.920232E-4, 4.461151E-4, 1.144283E-4, 
        5.926937E+0, 1.466600E-3, 1.947206E-3, 1.763130E-3, 4.720412E-4, 
        4.152104E+00, 2.938930E-3, 3.850049E-3, 3.485422E-3, 9.678104E-4, 
        3.696742E+00, 3.674192E-3, 4.787192E-3, 4.334504E-3, 1.221669E-3, 
        1.730869E+00, 1.448399E-2, 1.787203E-2, 1.631358E-2, 5.320333E-3, 
        1.127490E+00, 3.426555E-2, 3.290205E-2, 3.048440E-2, 1.125278E-2, 
        9.694371E-1, 5.096100E-2, 3.952124E-2, 3.689802E-2, 1.431383E-2, 
        2.870226E-1, 2.015616E-1, 9.850062E-2, 9.105633E-2, 5.377092E-2, 
        1.126476E-1, 1.323797E-1, 1.095493E-1, 9.780977E-2, 7.193724E-2, 
        7.712071E-2, 1.000004E-1, 9.660049E-2, 8.905284E-2, 7.049105E-2 },
        sabCorrect {1.193441E1, 3.645467E-4, 4.918339E-4, 4.459435E-4, 
        1.143843E-4, 5.917825, 1.464345E-3, 1.944212E-3, 1.760420E-3, 4.713155E-4, 
        4.139347, 2.929900E-3, 3.838220E-3, 3.474713E-3, 9.648392E-4, 3.682550, 
        3.660087E-3, 4.768813E-3, 4.317864E-3, 1.217011E-3, 1.704442, 1.426286E-2, 
        1.759921E-2, 1.606475E-2, 5.240240E-3, 1.093324, 3.322725E-2, 3.190517E-2, 
        2.956113E-2, 1.091342E-2, 9.328565E-1, 4.903809E-2, 3.803011E-2, 
        3.550622E-2, 1.377530E-2, 2.460920E-1, 1.728181E-1, 8.445423E-2, 
        7.807188E-2, 4.610450E-2, 8.281073E-2, 9.731646E-2, 8.053320E-2, 
        7.190342E-2, 5.288448E-2, 5.249625E-2, 6.807050E-2, 6.575628E-2, 
        6.061882E-2, 4.798440E-2};

      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, temp, oscEnergiesWeights, effectiveTemp, sab );

      THEN( "scattering law matrix is correctly changed" ){
        REQUIRE(ranges::equal(sab,sabCorrect,equal));
      } // THEN
    } // WHEN

    WHEN( "alpha and beta both span a huge range" ){
      std::vector<double> 
        alpha = {0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1.0, 5.0, 10, 50, 100},
        beta  = {0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1.0, 5.0, 10, 50, 100},
        sab {3.783747E+01, 3.389591E+01, 2.424696E+01, 5.632040E-04, 
        5.944319E-05, 2.931679E-05, 3.604359E-05, 1.302505E-05, 2.852185E-10, 
        0.000000E+00, 0.000000E+00, 1.690484E+01, 1.652874E+01, 1.552734E+01, 
        1.827156E+00, 2.486844E-03, 1.482647E-04, 1.794480E-04, 6.472644E-05, 
        7.025305E-09, 0.000000E+00, 0.000000E+00, 1.193902E+01, 1.182827E+01, 
        1.145869E+01, 3.965110E+00, 1.397768E-01, 2.990279E-04, 3.576237E-04, 
        1.289955E-04, 2.781304E-08, 0.000000E+00, 0.000000E+00, 5.288890E+00, 
        5.277051E+00, 5.265238E+00, 4.322478E+00, 2.256523E+00, 1.546854E-03, 
        1.755870E-03, 6.391400E-04, 6.727212E-07, 0.000000E+00, 0.000000E+00, 
        3.696938E+00, 3.693800E+00, 3.690662E+00, 3.379998E+00, 2.474803E+00, 
        3.230450E-03, 3.452499E-03, 1.274852E-03, 2.958873E-06, 5.010057E-30, 
        0.000000E+00, 1.515354E+00, 1.516883E+00, 1.518416E+00, 1.515707E+00, 
        1.457246E+00, 2.191122E-01, 1.593964E-02, 6.432143E-03, 8.011742E-05, 
        1.816200E-23, 0.000000E+00, 9.667891E-01, 9.682187E-01, 9.696371E-01, 
        9.799659E-01, 9.716895E-01, 4.202779E-01, 4.580241E-02, 1.304702E-02, 
        3.455963E-04, 1.457970E-20, 0.000000E+00, 2.128595E-01, 2.133201E-01, 
        2.137239E-01, 2.171325E-01, 2.218032E-01, 2.275933E-01, 1.776032E-01, 
        5.925186E-02, 1.045667E-02, 2.044901E-13, 0.000000E+00, 7.265244E-02, 
        7.282368E-02, 7.299636E-02, 7.441957E-02, 7.611087E-02, 8.707374E-02, 
        9.225780E-02, 7.963887E-02, 3.473308E-02, 3.651351E-10, 1.467265E-24, 
        4.117754E-04, 4.127892E-04, 4.138054E-04, 4.220265E-04, 4.325349E-04, 
        5.257059E-04, 6.655505E-04, 3.201093E-03, 1.151002E-02, 1.287332E-03, 
        6.291107E-10, 1.706488E-06, 1.710763E-06, 1.715049E-06, 1.749688E-06, 
        1.793778E-06, 2.187137E-06, 2.789660E-06, 1.680559E-05, 1.049101E-04, 
        1.531873E-02, 6.550072E-05}, 
        sabCorrect {3.783601E+01, 3.389461E+01, 2.424603E+01, 5.631824E-04, 
        5.944090E-05, 2.931567E-05, 3.604220E-05, 1.302457E-05, 8.703740E-10, 
        0.000000E+00, 0.000000E+00, 1.690159E+01, 1.652556E+01, 1.552435E+01, 
        1.826804E+00, 2.486366E-03, 1.482362E-04, 1.794134E-04, 6.471452E-05, 
        2.158185E-08, 0.000000E+00, 0.000000E+00, 1.193442E+01, 1.182372E+01, 
        1.145428E+01, 3.963585E+00, 1.397231E-01, 2.989129E-04, 3.574861E-04, 
        1.289481E-04, 8.581651E-08, 0.000000E+00, 0.000000E+00, 5.278728E+00, 
        5.266912E+00, 5.255121E+00, 4.314173E+00, 2.252187E+00, 1.543882E-03, 
        1.752496E-03, 6.379641E-04, 2.096558E-06, 1.025453E-18, 0.000000E+00, 
        3.682745E+00, 3.679619E+00, 3.676493E+00, 3.367022E+00, 2.465302E+00, 
        3.218048E-03, 3.439245E-03, 1.270163E-03, 8.560432E-06, 1.821032E-13, 
        0.000000E+00, 1.486489E+00, 1.487989E+00, 1.489492E+00, 1.486835E+00, 
        1.429487E+00, 2.149384E-01, 1.563602E-02, 6.314523E-03, 2.089396E-04, 
        2.440142E-10, 0.000000E+00, 9.303083E-01, 9.316840E-01, 9.330489E-01, 
        9.429879E-01, 9.350238E-01, 4.044191E-01, 4.407416E-02, 1.257785E-02, 
        1.008290E-03, 5.611201E-08, 0.000000E+00, 1.756196E-01, 1.759995E-01, 
        1.763327E-01, 1.791450E-01, 1.829986E-01, 1.877760E-01, 1.465330E-01, 
        4.930239E-02, 2.029639E-02, 3.243924E-05, 2.539787E-12, 4.945904E-02, 
        4.957562E-02, 4.969317E-02, 5.066204E-02, 5.181344E-02, 5.927721E-02, 
        6.280896E-02, 5.479071E-02, 3.622037E-02, 3.282077E-04, 7.837533E-08, 
        6.087852E-05, 6.102842E-05, 6.117870E-05, 6.239442E-05, 6.394848E-05, 
        7.773404E-05, 9.845356E-05, 4.792199E-04, 1.831175E-03, 3.620962E-03, 
        3.530294E-04, 3.795997E-08, 3.805508E-08, 3.815042E-08, 3.892100E-08, 
        3.990205E-08, 4.865977E-08, 6.209764E-08, 3.764874E-07, 2.441230E-06, 
        8.235741E-04, 3.394048E-03};

      discreteOscillators( lambda_s, twt, tbeta, alpha, beta, temp, oscEnergiesWeights, effectiveTemp, sab );

      THEN( "scattering law matrix is correctly changed" ){
        REQUIRE(ranges::equal(sab,sabCorrect,equal));
      } // THEN


    } // WHEN
  } // GIVEN

} // TEST CASE
